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ABSTRACT 

During a galaxy merger, the supermassive black hole (SMBH) in each galaxy is thought to sink to 
the center of the potential and form a supermassive black hole binary; this binary can eject stars via 3- 
body scattering, bringing the SMBHs ever closer. In a static spherical galaxy model, the binary stalls 
at a separation of about a parsec after ejecting all the stars in its loss cone - this is the well-known 
final parsec problem. Earlier work has shown that the centrophilic orbits in triaxial galaxy models are 
key in refilling the loss cone at a high enough rate to prevent the black holes from stalling. However, 
the evolution of binary SMBHs has never been explored in axisymmetric galaxies, so it is not clear if 
the final parsec problem persists in these systems. Here we use a suite of direct A^-body simulations 
to follow SMBH binary evolution in galaxy models with a range of ellipticity. For the first time, we 
show that mere axisymmetry can solve the final parsec problem; we find the the SMBH evolution is 
independent of N for an axis ratio oi c/a — 0.8, and that the SMBH binary separation reaches the 
gravitational radiation regime for c/a = 0.75. 

Subject headings: Stellar dynamics - black hole physics - Galaxies: kinematics and dynamics - Galaxy: 
center. 



1. INTRODUCTION 

Supermassive Black Holes (SMBHs), with masses 
from 10^ — 10^°Mq are thought to dwell at the cen- 
ters of most galaxies (e .g. iKormendv fc Richstond [19951 : 
iFerrarese fc FordI 120051 ). and since we know that galaxy 
mergers are commonplace, the formation of a binary 
SMBH is almost inevitable at so me point during galaxy 
assembly ([Begelman et al.l 119801 ). There are cases in 
which there is a clear observ ational evidence for tw o 
widely separated SMBHs(e.g. iComerford et al.l l2009f ). 
as w ell as growing evidence f or true SMBH bina- 
ries (Komossa] 120061 : iBonI 120121 ). Upon coalescence, 
SMBH binaries would be the highest signal-to-noise ra- 
tio sources of gravitational waves, detectable from es- 
sentially any cosmological redshift, for future space- 
borne gravitational wave detectors such as the New 
Gravitational wave Observatory CNGO) (IHughesI 12001 



Harack & Cutler 2004; Danzm annet al.ll201lD 

The paradigm for SMBH binary evolution consists of 
three distinct phases ( Begelman et al.lll980D . First, the 
two SMBHs sink towards the galaxy center due to the 
dynamical friction, eventually forming a bound pair with 
semi-major axis a ^ rji, where is the binary's influence 
radius, typically is defined as the radius enclosing twice 
the mass of the binary in stars M{< r^) = 2M,, where 
M, is the mass of the binary. Dynamical friction contin- 
ues to drive the SMBHs closer until the system becomes 
a hard binary; here the separation is a ~ ah : (jQuinlanl 

IT99I r^[2ooa) 



Oh 



4^ 



(1) 



where is the binary's reduced mass and a is the ID ve- 
locity dispersion. In this second phase, slingshot ejection 
of stars is the dominant mechanism for taking energy and 



angular momentum away from SMBH binary. If 3-body 
scattering can shrink the binary orbit by a factor of or- 
der a hundred, then the SMBHs will be close enough to 
emit copious gravitational radiation. This gravitational 
wave emission in the third and final phase will drain or- 
bital energy from the binary and will drive the final co- 
alescence in roughly 100 Myr. The transition from the 
first to the second phase can be quite prompt, provided 
that the ga laxy mass ratio is not too small (g > 0.1) 
(jCallegari et al . 2011). In contrast, the subsequent tran- 
sition from the 3-body scattering to gravitational radia- 
tion regime could constitute a bottleneck for the binary 
coalescence, simply caused by a lack of stars which can 
interact with the binary. This is the so-called Final Par- 
sec Problem. 

Theorists commonly invoke gas to move these SMBHs 
through the final parsec. This is because the SMBHs can 
frict against the gas ver y efficiently, as long as the gas 
is no t forming stars fe.g lMaver et al.ir2007l : lEscala et al.l 
l2005f ). Unfortunately, relying on gas to solve the Final 
Parsec Problem will leave out the most massive SMBHs, 
since their host ga laxies are massive ellipticals and are 
therefore gas-poor (|Kannappanl [20Ql). 

There has been significant work to describe the Fi- 
nal Parsec Problem and to solve it with a purely stel- 
lar dynamical approach. For example, lYul (|2002D first 
pointed out, based on an analysis of an HST sam- 
ple of nearby ellip tical galaxies and spiral bulges from 
iFaber et al.l (|1997l ). that flattening and non-axisymmetry 
could bring SMBH binaries to coalescence. Later, sev- 
eral studies built self-consistent cuspy, triaxial mod- 
els with a single central SMBH and showed that 
models could be concocted with a signiflcant fraction 
of centrophilic orbits that would efficiently drive the 
hardening rate of a binary if it were present (e.g. 
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Merritt fc Poot3 120041: iHollev- Bockelmann et all 120021: 
HoUev-Bockelmann fc S igurdssonl l2006[ ). Also direct AC- 



TABLE 1 

Parameters of the SMBH binary study 



body simulations of SMBH binary evolution in iso- 
lated galaxy models showed that the Final Parsec Prob- 
lem could potentially be solved by purely stellar dy- 
namical m odels that develop a significant amount o f 
triaxiality (jBerczik et all 120061; iBerentzen et alj I2OO90 . 
Recent numerical studies showed th at galaxies which 
form via mergers are mildly triaxia l ()Khan et alJ 120111 : 
iPreto et alJ 120111 : iKhan et all l2012f ). Hardening rates 
of SMBH binaries resulting from equal-mass mergers 
are substantially higher than those found in spherical 
nuclei, and are independent of the total number N of 
stars. Presumably, the triaxial shapes of merger rem- 
nants reported in these studies result in a large popula- 
tion of stars on regular centrophilic orbits deep within 
the galaxy potential, such as pyrami d or box orbits with 
form a lly zero angular momentum IG e rhard fc Binnev 
19851: iValluri fc MerrittI 119981: IHollev -Bockelman n et al 



2001f ). Outside the influence radius of the SMBH, many 
of the centrophilic orbits are chaotic in a triaxial po- 
tential . th ough some regular resonant boxlet or bits re- 
main (| Po on & Merritt 2001; HoUcv-Bockelman n et alJ 
[2001 . 

In contrast to the wealth of work on solving the final 
parsec problem in triaxial galaxies, comparatively little 
attention has been paid to exploring binary coalescence 
in axisymmetric systems. Recent work on the capture 
of stars from the loss cone in axisymmetric nuclei looks 
promising for the final parsec problem, in that the stel- 
lar ca pture rate is several tinies higher in flattened sys- 
tems (jMagorrian fc Tremaind Il999t iVasiliev fc Merritt 
|2013[ ). In an axisymmetric galaxy, the centrophilic orbit 
family includes saucer or cone orbits (jSridhar fc Toumal 
[1999) ■ and these can, in principle increase 3-body scat- 
tering rates. However, no detailed N-body simulations 
have been done to test this. In this paper, we explore 
SMBH coalescence in axisymmetric galaxies using high 
resolution N-body simulations of equilibrium axisymmet- 
ric galaxy models with a range of intrinsic flattening. 

2. INITIAL CONDITIONS AND NUMERICAL METHODS 

2.1. The host galaxies and their SMBHs 

We generated our equilibrium axisymmetric sys- 
tem s from initially spheri cal SMBH-embedded rj mod- 
els (jTremaine et al.l 119941 ) by applying a slow, steady 
velocity drag on the particles in the z direction (see 
IHollev-Bockelmann et all 120011 for details). This 
adiabatic squeezing technique preserves the phase 
space distribution function of the system, so al- 
though the orbit content changes, the density pro- 
file rem ains unchanged as the galax y shape slowly 
evolves (jHollev-Bockelmann et alll2002f ). The suite we 
present here has an inner density slope, 7, of 1.0. Our 
model has a mass of 1.0 in system units, the central 
SMBH has a mass of 0.005, and the half mass radius is 
2.41; scaling the SMBH mass to the Milky Way implies 
that we are simulating the inner few hundred parsecs of 
the galaxy. A second equal mass SMBH is introduced 
into the outskirts of the equilibrium model at a distance 
of 0.5 in model units with 70 % of the circular velocity. 
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Note. — Column 1: Galaxy model. Column 2: Number of 
particles. Column 3: Central density slope 7. Column 4; Axes 
ratio. Column 5: SMBH mass ratio. 

The A^-body integrations are carried out using 
the 0GRAPEI3 (jHarfst et al.l I2007D . a parallel, direct- 
summation A^-body code that uses Graphic Processing 
Units (GPU) cards to accelerate the computation of pair- 
wise gravitational forces between all particles. 0GRAPE 
code uses a 4th-order Hermite integration scheme to 
advance all particles with hierarchical individual block 
time-steps, together with the parallel usage of GPU cards 
to calculate acceleration a and its first time derivative a 
between all p articles. For the forc e calculation we use 
the SAPPORO (jGaburov et al.l 120091 ) library which emu- 
lates standard GRAPE-6 library calls on the GPU. 

As (/)GRAPE does n ot include the regularization 
(jMikkola fc AarsethI Il998f ) of close encounters or bina- 
ries, we use a gravitational Plummer-softening for force 
calculations. The softening length for star-star encoun- 
ters is chosen to be smaller (i.e. 10~^ in our model units) 
than the minimum separation reached by SMBH binaries 
in our study. The softening between the SMBH particles 
is set equal to zero. 

The suite of A^-body experiments were carried out on 
a high-performance GPU computing cluster ACCRE em- 
ploying 192 GPUs at Vanderbilt University, Nashville, 
TN. 

For a detailed discussion on the A^-body methods and 
choice of parameters, we refer the reader to IKhan et al.l 
((20II . 

3. SMBH BINARY EVOLUTION IN FLATTENED GALAXY 
MODELS 

We created models with axis ratios c/a = 
[1.0,0.9,0.8,0.75] (see table [T] ) in order to study the 
SMBH binary evolution as a function of flatness. 

Figure [1] describes the comparison between the evo- 
lution of the SMBH binary in spherical and flattened 
galaxy models (c/a = 0.9). In the first phase, which 
can be seen before T = 10, the separation between the 
two SMBHs shrink due to dynamical friction. The evo- 
lution of SMBHs in this phase is very similar for both 
the spherical and flattened galaxy models. However the 



2.2. Numerical Methods and Hardware 
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Fig. 1. — Evolution of the separation between each SMBH 
(top), the semimajor axis (middle) and eccentricity (bottom) of the 
SMBH binary in Af-body integrations for spherical and flattened 
models with c/a = 0.9 (see table[TJ. Time and R are measured in 
model units. 
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Fig. 2. — Average hardening rates of SMBH binary in spherical 
and flattened galaxy models in Table [T] The average is measured 
between T = 40 and T = 70 in model units. 

in both of these models. Note, though, that s is roughly 
four times higher for the SMBH binary evolution in flat- 
tened system (c/a — 0.9) compared to a spherical galaxy 
model. In flattened galaxy models, then, there appears 
to be an additional supply of stars that interact with the 
massive binary. 

The SMBH binary evolution for non-rotating models 
with c/a ~ 0.8 are shown in figure [H On the face of it, 
the evolution of SMBH binary in flattened models with 
c/a = 0.8 is very similar to c/a = 0.9. However, there 
is a striking difference for runs with N greater than 500k 
in the sense that SMBH binary evolution becomes inde- 
pendent of N. This implies that the system is no longer 
dominated by artificial 2-body relaxation, and we are 
uncovering the accurate SMBH binary evolution. Again, 
we calculated SMBH binary hardening rates in interval 
between T = 30 and T = 60. Note that for N > 500k, 
s = 22 is 8 times higher than the spherical galaxy model 
at the highest particle resolution. 

Figure [4] features the SMBH binary evolution for the 
flattest model we explored so far, c — 0.75, and we con- 
firm the both N-independence and the rapid hardening 
rate seen in the c = 0.8 runs. As is perhaps expected, 
the hardening rate is slightly enhanced compared to the 
c = 0.8 model and the SMBHs reach an even smaller 
separation. 



subsequent evolution, governed by three body encoun- 
ters, happens at much faster rate in flat galaxy models, 
as can be seen from the evolution of both the separation 
and the semimajor axis of the binary SMBH (the top 
and middle panels). The SMBH binary forms with an 
eccentricity, e, much less than 0.1 and it remains small 
during the binary evolution. Although the semi-major 
axis of SMBH binary evolves faster than that of spher- 
ical galaxy model, it still shows clear dependence on N 
(see middle panel) - a sure sign that the binary evolution 
in this simulation is plagued by the numerical effects of 
2-body relaxation. 

The binary hardening rates s are calculated by fitting 
straight lines to a~^{t) in the late phase of binary evo- 
lution from T = 30 — 60. The value of s decreases with 
increasing iV, and again this suggests that numerical 2- 
body relaxation is important for SMBH binary evolution 



3.1. Estimates for Relativistic Regime 

At the end of our simulations, the SMBH binary 
separation is still far larger than what is required for 
gravitational wave emission to dominate and guide the 
binaries on their way to coalescen ce. Earlier studies 
(jBerentzen et all 120091 : iKhan et all (2012) have shown 
that the estimated coalescence time obtained using con- 
stant hardening rat e s in t he ste llar dynamical regime 
and the formula of iPetersI ()1964[ ) for hardening in the 
gravitational wave dominated regime agree remarkably 
well with Tcoai obtained from simulations that follow the 
binary evolution until coalescence using post-Newtonian 
[VN) terms in the equation of motion of the SMBH bi- 
nary. 

So, we will use the lPeteri ()1964[ ) formalism to estimate 
the evolution o f the SMBH binar y beyond the Nbody 
regime (see also lKhan eralll2012af) : 
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Fig. 3. — Evolution of the binary SMBH separation (top), semi- 
major axis (middle) and eccentricity (bottom) of the SMBH binary 
in A'^-body integrations for spherical and flattened models with 
c/a = 0.8 (see table [T)l . Time and K are measured in model units. 
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304' 



(4) 



da 



NB 



/da\ 
\dt / GW 



+ (2) 
\dt / GW 



The orbit-averaged expressions — including the lowest 
order 2.5 VM dissipative terms — for the rates of change 
of a binary's semi-major axis, and ec centricity due to 
gravitational wave emission are given bv lPetersl(fl96l : 



da\ 

dt I GW 



I de\ 

\ dt I GW 



64 G'3M.iM.2(M.i + Af.2) 

h a3c5(l -e2)7/2 

73 2 37 4 

H H 

24 96 

304 G'3M.iM.2(M.i + M,2) 

a4c5(l„e2)5/2 ^ 



(3) 



In order to compute the evolution in the gravitational 
radiation regime, we must adopt some physical length 
and mass. We use the observed mass of the SMBH and 
its influence radius for NGC4486A and compare it to the 
mass and influence radius of our model. One unit of 
length and mass in our model units equals 0.3 kpc and 
2.6 X 10^ solar masses, respectively. The value of speed 
of light is c = 1550 and one time unit is equal to 1.4 Myr. 

We next present an estimate of the coalescence time 
for the run Flat75a, which is a typical representation 
of the c/a = 0.75 runs. At T = 70 in model units, the 
semi-major axis a is 5 x 10~^ and the eccentricity e is 0.1, 
which makes the estimated coalescence 2.4 Gyr. This co- 
alescence is a few times longer than the coalescence times 
for radial galaxy mergers in , Khan et al.. (,2012i) . but it is 
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still well within a Hubble time. We speculate that the 
main reason behind this difference is that the eccentric- 
ity here is quite low. Indeed, we are currently exploring 
the eccentricity evolution in axisymmetric galaxy mod- 
els to determine if low eccentricity is a generic feature of 
flattened equilibrium models. 

4. SUMMARY & CONCLUSIONS 

We present the first results of SMBH binary evolu- 
tion in a fully self-consistent equilibrium axisymmetric 
galaxy model. We find, in general, that SMBH bina- 
ries evolve faster and reach smaller separations for even 
mildly flattened galaxy models. For c/a = 0.8, we begin 
to see hints of TV-independence for > 500fc, and here 
the hardening rate is 25 times faster than in a spherical 
model. For our flattest model c/a — 0.75, the system 
evolves to a gravitational radiation regime regardless of 
the particle resolution; for sufhciently flattened systems, 
we show that the final parsec problem is not a prob- 
lem. These models may best pertain to pseudobulges 
and SOs; a recent example is the lenticular galaxy NGC 
1277, whic h seems to host a new clas s of ultramassive 
black hole (jvan den Bosch et al.l 12012} ). Note, though, 
that the flattening in our models is very modest com- 
pared to these syste ms, where the axis ratios are often 
c/a ~ 0.6 or smaller (jKormendv fc Kennicuttl[2004D . and 
wc have also yet to include rotation. Both of these effects 
should act to increase the hardening rate. 

Contrary to previous work that consistently predict 
highly eccentric SMBH binaries as they pass to the grav- 
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itational radiation regime (IKhan et al.ll201lL[20Tlal ). we 
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centricity, to unravel the eccentricity evolution of these 
SMBH binaries. 

In this pilot study, all of our flattened galaxy models 
had the same density profile, central black hole mass, 
black hole mass ratio, and initial black hole orbital ec- 
centricity. Since the SMBH binary evolution will likely 
depend on the internal structure of the host galaxy and 
on details of the black hole orbit, we are launching a 
systematic study of SMBH binary evolution in axisym- 
metric galaxy models of various cusp slopes, mass ratios, 
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This work was conducted at the Advanced Comput- 
ing Center for Research and Education at Vanderbilt 
University, Nashville, TN. K. H.-B. acknowledges sup- 
port from NSF CAREER award AST-0847696, as weU 
as the support from the Aspen Center for Physics. K.H.- 
B. also wishes to thank Daniel Sissom for assistance in 
this project. 



Khan, F. M., Just, A., & Merritt, D. 2011, ApJ, 732, 89 
Khan, F. M., Preto, M., Berczik, P., Berentzen, I., Just, A., & 

Spurzem, R. 2012, ApJ, 749, 147 
Khan, F. M., Berentzen, I., Berczik, P., Just, A., Mayer, L., 

Nitadori, K. & Callegari, S., 2012, ApJ, 756, 30 
Komossa, S. 2006, Memorie della Societ Astronomica Italiana, 77, 

733 

Kormendy, J., & Kennicutt, R. C, Jr. 2004, ARA&A, 42, 603 
Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581 
Magorrian, J., & Tremaine, S. 1999, MNRAS, 309, 447 
Mayer, L., Kazantzidis, S., Madau, P., et al. 2007, Science, 316, 
1874 

Merritt, D., & Poon, M. Y. 2004, ApJ, 606, 788 

Merritt, D., & Vasiliev, E. 2011, ApJ, 726, 61 

Mikkola, S., Aarseth, S. 1998, New Astronomy, 3, 309 

Peters, P. C. 1964, Phys. Rev. B, 136, 1224 

Poon, M. Y., & Merritt, D. 2001, ApJ, 549, 192 

Preto, M., Berentzen, I. Berczik, P., Merritt, D. & Spurzem, R. 

2009, Journal of Physics: Conference Series, Volume 154, Issue 

1, pp. 012049 

Preto, M., Berentzen, I. Berczik, P., & Spurzem, R. 2011, ApJ, 
732, L26 

Quinlan, G. D. 1996, New Astron., 1, 35 
Sridhar, S., & Touma, J. 1999, MNRAS, 303, 483 
Tremaine, S., Richstone, D. O., Byun, Y.-L, et al. 1994, AJ, 107, 
634 

Valluri, M., & Merritt, D. 1998, ApJ, 506, 686 

van den Bosch, R. C. E., Gebhardt, K., Giiltekin, K., et al. 2012, 

Nature, 491, 729 

Vasiliev, E., & Merritt, D. 2013, arXiv:1301.3150l 
Yu, Q. 2002, MNRAS, 331, 935 



URL (cited on August 2011) pittp: //sci.esa.int/science-e/www/object/index.cfm?fobjectid=48364| 



